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Abstract 

This paper describes a general-purpose programming technique, called the Simulation of 
Simplicity, which can be used to cope with degenerate input data for geometric algorithms. 
It relieves the programmer from the task to provide a consistent treatment for every single 
special case that can occur. The programs that use the technique tend to be considerably 
smaller and more robust than those that do not use it. We believe that this technique will 
become a standard tool in writing geometric software. 
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1 Introduction 

This paper introduces a general technique that can be used to cope with degenerate cases encoun- 
tered by computer programs. Consider, for example, a program that sorts an array of integers 
using a comparison as a primitive operation. A special, or degenerate, case occurs when the pro- 
gram attempts to decide which one of two equal numbers is smaller than the other. A typical 
way to resolve this tie is to pretend that the number with smaller index is smaller (assuming 
the integers are indexed, e.g., by their positions in an array). Or think of Kruskal's algorithm 
for constructing a minimum spanning tree of a weighted graph (see [AHU74]). At each step it 
chooses the shortest edge that can be added to the current collection of edges without creating a 
cycle. If this edge is not unique, then any one of the candidate edges is taken. The thus generated 
minimum spanning tree is therefore not unique unless we specify deterministic rules to break ties. 

In both problems, sorting and constructing minimum spanning trees, the special cases are easily 
dealt with, partly because the ties can be broken arbitrarily without creating inconsistencies. 
The situation is usually far more complicated for geometric problems. Consider for example the 
following seemingly straightforward algorithm for the point-in-polygon problem which is sometimes 
called the Parity Algorithm. 

• Let r be the horizontal half-line whose left endpoint is the test point. 

• Count the number of intersections between r and the edges of the polygon. If that number 
is odd, then the test point lies within the polygon, and if the number is even, then it lies 
outside the polygon. 

As pointed out in [Fo85], it is not a trivial matter to implement this algorithm, even if we assume 
that the test point does not lie on the boundary of the polygon. There are only two nondegenerate 
cases: Either the intersection between r and an edge e is empty or r crosses e (see Figure l-I, (a) 
and (b)). There are, however, four degenerate cases (as illustrated in Figure l-I, (c) through (f)) 
that have to be taken into account. 



A correct answer is obtained if cases (c) and (e) are counted as one crossing and cases (d) and (f ) 
are not counted at all. If we write the code for the above algorithm, we realize that a substantial 
amount of the effort is required to cover the four degenerated cases. Observe also that there are 
several seemingly plausible ways to treat the degenerate cases and that some of them lead to 
incorrect algorithms. We appeal to the imagination of the reader to envision the bizarre structure 
of degenerate cases one encounters in generalizing the point-in-polygon problem to three or higher 
dimensions. Another problem with a set of degenerate cases that is considerably richer than the 
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one of the point-in-polygon problem is obtained if one intersects a polygon with a geometric object 
that is more complicated than a half-line. 

When it comes to implementing geometric algorithms, degenerate cases are very costly, in partic- 
ular if there are many such cases that have to be distinguished. This is caused by the positive 
correlation between the number of degenerate cases and a variety of factors that contribute to the 
overall cost of a piece of software. These factors include the length of the program, which, for itself, 
correlates positively with the amount of time required to write it, to debug it, and to maintain it. 
Of course, the degree of robustness of the program decreases with increasing complication. The 
correctness of a program relies on the consistent treatment of all different cases. In this context, 
it is worthwhile to mention that more efficient algorithms tend to be more complicated and also 
more sensible to slight inconsistencies in treating degenerate cases. 

This paper presents a general technique, called Simulation of Simplicity (SoS), that can be used to 
cope with the problems mentioned above. Intuitively, it simulates a conceptual perturbation of the 
input data that eliminates all degeneracies. We hasten to mention that the perturbation is never 
ever computed — it is assumed to be arbitrarily small, although not vanishing, which is enough to 
simulate the nondegenerate topology. Another interpretation of the technique views it as a general 
way to break ties consistently. The tie-breaking part of the code appears in the lowest level of the 
algorithm, namely, in the procedures that implement the needed primitive operations. Different 
techniques following the same main approach have recently been suggested in [Ya87, Ya88]. A 
large part of this paper is devoted to demonstrating that the overhead in time caused by the use 
of the more elaborate primitive procedures required by SoS is negligible. 

The outline of this paper is as follows: Section 2 presents the general idea of the technique 
and works out some guidelines needed to implement it effectively. Section 3 considers a class of 
problems for finite point sets that can be solved using a common set of geometric primitives. It 
also discusses how the perturbation infiuences the geometric primitives. Section 4 demonstrates 
efficient implementations of the primitive operations. In Section 5 we show that the geometric 
primitives introduced for point set problems can be used to solve a variety of other problems defined 
for polygons, hyperplanes, circles, spheres, and other geometric objects. Finally, in Section 6, we 
discuss the perturbation technique and its limitations. 

2 SoS — the General Idea 

Degeneracies occur with probability zero if we draw a finite number of geometric objects, each 
represented by a finite set of numbers from the (infinite) set of all such objects, provided there is 
no bound on the precision of the numbers used. In real-life computing this is not the case; that 
is, there is only a finite set of available numbers and thus a bound on the precision that can be 
achieved. As a consequence, we are doomed to work with degenerate data. On the other hand, 
even infinite precision does not guarantee the nonexistence of degeneracies. This section gives the 
general outline of a technique called the Simulation of Simplicity (SoS) — we use "simple" as a 
synonym for nondegenerate — which allows us to neglect degeneracies when we write programs. 
A similar but less elaborate method has been used to solve degenerate linear programs. This leads 
to the implementation of the simplex algorithm referred to as the "lexicographical method" (see 
[Ch52], [DOW55], [Da63], or [Ch83] for details). In computational geometry, this technique has 
been used in a couple of papers, including [Ed86] and [EW86], to avoid the otherwise necessary 
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discussion of degenerate cases. This paper presents the theoretical foundations of SoS as well as 
details of its implementation. 

The basic idea of SoS is to perturb the given objects slightly which amounts to changing the 
numbers that represent the objects; these numbers will be called the coordinates or the parameters 
of the objects. It is important that the perturbation is small enough so that it does not change 
the nondegenerate position of objects relative to each other. Coming up with such a perturbation 
is rather difficult and may require much higher precision than used for the original set of objects. 
For this reason, we perform the perturbation only symbolically by replacing each coordinate by a 
polynomial in e. The polynomials will be chosen in such a way that the perturbed set goes towards 
the original set as e goes to zero. We will see that it is not important to know the exact value 
of e to perform the simulation; rather, it is sufficient to assume that e is positive and sufficiently 
small. Thus, it will be possible to use e as an indeterminant and to handle primitive operations 
symbolically. 

The future user of SoS will neither have to be concerned with the role that e plays in the perturba- 
tion nor with the symbolic manipulation of polynomials. We may think of SoS as a package that 
provides the primitive operations needed for a certain computation. Ideally, the inside of these 
operations is hidden from the user who communicates with them like with an oracle. It turns out 
that a large number of geometric problems can be solved using a surprisingly small number of 
primitives. Some of these primitives will be discussed in the following three sections. This section 
continues to develop the general ideas on which SoS is based. 

One of the goals of SoS is to perturb a set of objects such that all degeneracies disappear. A 
degeneracy is something that is not defined in general; its definition depends on the problem at 
hand. More specifically, it depends on the primitive operations used to solve the problem. For 
example, a primitive operation in the point-in-polygon algorithm described in the introduction 
tests the intersection of a horizontal half-line and a line segment. A degeneracy occurs if the 
half-line contains one or both endpoints of the line segment. A set of objects is now called simple^ 
or nondegenerate^ or in general position^ if it does not contain any degeneracy. We thus define 
"simplicity" relative to the primitives used to solve a problem. 

This paper considers only topological primitives^ that is, operations that test some given input and 
classify it as one of a constant number of possible cases. This is in contrast to operations that 
compute new objects such as the intersection of a half-line and a line segment. In most programs, 
such an object serves only as an intermediate result anyway; but an intermediate result can as 
well be represented implicitly as a collection of pointers and a tag that tells us in what sense the 
objects identified by the pointers determine the (implicit) result. To simplify our discussion even 
further, we restrict our attention to primitives with three possible outcomes which we represent 
by -|-I, 0, and —I, where indicates a degeneracy and -|-I and —I distinguish between the two 
nondegenerate cases. Tests that distinguish between more than two nondegenerate cases can be 
obtained by combining several ternary tests. 

If we think of a primitive operation as a function / that maps a high-dimensional point (whose 
coordinates describe the input objects) to -|-I, 0, or —I, then /~^(0) represents the set of degenerate 
inputs. One requirement for this set is that its measure in this high-dimensional space is zero — 
otherwise, it is unreasonable to call its points degenerate. A set of n objects, given by d parameters 
each, can be thought of as a point in nd dimensions. If / takes k < n objects as input, then /~^(0) 
is a surface of measure zero in A;J-dimensional space. This surface defines another zero-measure 
surface in nd dimensions which is obtained by embedding /~^(0) in the A;J-dimensional subspace 
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defined by the k objects and extending it orthogonal to this subspace along the other coordinate 
axes. Other combinations of k objects provide additional zero-measure surfaces that, altogether, 
decompose the nJ-dimensional space into faces of various dimensions. A cell is an nJ-dimensional 
face of this decomposition, and all points of a cell correspond to nondegenerate sets of objects. 
A degenerate set corresponds to a point x in the union of the surfaces, denoted by S. Since S 
has measure zero, every nonempty open ball around this point contains a point y of some cell. 
Moving X to y corresponds now to perturbing the set of objects that x corresponds to such that all 
degeneracies disappear. This shows that a perturbation to a nondegenerate set is always possible 
even if the amount of perturbation is severely limited. Recall that another requirement for the 
perturbation is that it does not change any nondegenerate subconfiguration. This means that we 
should not move x across a surface it did not belong to initially. This can always be guaranteed if 
we choose the open ball small enough that it does not intersect any surface that does not contain 
the initial position of x. 

To follow the forthcoming reasoning it is not necessary for the reader to understand the topology 
of the nJ-dimensional space as indicated in the above paragraph. Nevertheless, this view of the 
problem sheds some light on the nature of degeneracy. It also explains why there is always a small 
enough perturbation that removes all degeneracies. Below, we discuss such perturbations more 
specifically and address a few questions concerning the efficient implementation of SoS. 

Simplicity is simulated by applying a particular perturbation to a set P = {po,Pi, • • • tPu-i} of n 
geometric objects 

Pz = (tTj,!, ^^^,2, • • • , T^z,d), < I < u - 1, 

each specified by d parameters. It will be important that each object has a unique index between 
and n — 1. The objects are in arbitrary, and therefore not necessarily in simple, position. The 
perturbation of P is realized by replacing each parameter by a polynomial in e. We define 

P(e) = {Pi{e) = (7r^,i(e), 7r^,2{e), • • • , t^zA^)) I < « < n - 1}, 

where 

7rij(e) = TTij + e{i,i) for 0<z<n — 1, l<i<c^, 

and e{i,i) a polynomial in e that goes to zero when e goes to zero. We will refer to the new 
parameters 7rij(e), the new objects Pi(e), and the new set P{e) as the e-expansions of the original 
parameters tTj-j, the original objects p^, and the original set P, respectively. The choice of the 
polynomials e(z,j) will be guided by three requirements SoS has to meet. 

(a) P(e) must be simple if e > is sufficiently small. 

(b) P(e) must retain all nondegenerate properties of the original set P. 

(c) The computational overhead caused by simulating P(e) should be negligible. 

As mentioned before, condition (b) is automatically met if e is small enough. To satisfy (a), it is 
sufficient to choose the e(z,j) such that there is no nonempty open interval / with the property 
that P(e) is not simple if e G /. Think of P as a point x in nd dimensions and let x{e) be the 
point that corresponds to P(e). The points a;(e), e > 0, form a one-dimensional curve C in nd. 
dimensions. Thus, (a) is satisfied if C fl iS is a discrete set of points. (Recall that S represents 
all points in nd. dimensions that correspond to degenerate sets P.) In this topological setting, the 
phrase "e sufficiently small" gets a specific meaning: If Eq > is the smallest value of e such that 
x{eo) G <S, then e is sufficiently small if and only if < e < Eq. It is less clear how condition (c) 
infiuences the choice of the e(z, j). Below, we formulate a criterion for the polynomials e(z, j) that 
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leads to an efficient implementation of SoS. However, we do not claim that other choices of the 
e(z,j) cannot lead to efficient implementations too. 

Recall that a primitive operation is a function / that maps a set Q of k objects to +1, 0, or — 1. 
If the e-expansion is defined properly, then f[Q[e)) G { + 1, —1} provided e > is small enough. 
In general, f[Q[e)) will be the sign of a fairly complicated function in e. (Since / is now a binary 
function we can identify { + 1, —1} with {true, false} and express it as a predicate. We will follow 
this practice in the following sections of this paper.) One way to allow for an efficient evaluation 
of f[Q[e)) is to choose the e{i,i) in different orders of magnitude such that two expressions, each 
consisting of several factors of the form e(z,j), can be compared solely on the basis of the index 
pairs involved. When we evaluate /((^(e)), we can sort its terms in order of decreasing 

significance which can be done by comparing sets of index pairs. The most significant term will 
be a term without any e-factor; it will be equal to f{Q). The first term with a nonzero coefficient 
decides the sign of the function. If Q is nondegenerate to begin with, then f[Q[e)) = f{Q)^ and 
no other term has to be determined. In Sections 3 and 4, we will see that such a choice of the 
e(z,j) allows us to determine the sign of a fairly complicated polynomial in only a few steps. 

Note that SoS requires us to tell when Q is degenerate, which means that we need to be able to 
decide whether or not f{Q) = 0. This is not possible with the kind of fioating-point arithmetic 
that it is usually provided by current computers. Instead, we need to use exact arithmetic and, 
thus, occasionally long integers. These admittedly somewhat expensive operations occur only 
inside the primitives and do not concern the user of SoS. Furthermore, the length of such long 
integers is bounded by a constant if kd^ the number of input parameters of /, is bounded by a 
constant. In most geometric algorithms, this constant is reasonably small. In Section 6 we report 
on our experience in implementing SoS and give an indication to what extent the use of long 
integer arithmetic slows down the computation. This point cannot be taken lightly because the 
long integer arithmetic is likely to occur in the innermost loop of any program that uses SoS and 
thus dictates the constant in front of the asymptotic running time. However, it is worthwhile to 
mention that the need for exact arithmetic is not a peculiar feature of SoS itself, but is necessary 
whenever we do exact computation rather than push our luck and hope for the cancellation of 
round-off errors. 

3 Finite Point Sets — a Case Study 

For a further discussion of SoS it is advantageous to apply it to certain geometric objects and 
certain primitive operations defined for these objects. We choose points in the J-dimensional 
Fuclidean space E'^ as the objects for the case study. Notice that this is actually no loss of 
generality since every object specified by d parameters can be interpreted as a point in E'^. The 
primitive operation that we will consider takes d -\- 1 points as input and decides on which side of 
the hyperplane spanned by the last d points the first point lies. As we will see in Section 5, this 
primitive operation has a wide range of applications. 

If a given finite point set is perturbed, as explained in Section 2, one can ignore all degeneracies 
and special cases. The price for this simulated simplicity is that the coordinates of the points are 
now symbolic expressions in e. Fven for a simple task, such as the comparison of two coordinates, 
we need a custom-made procedure that handles the e-expansions of the coordinates. Let tTj-j be 
the j-th coordinate of point pi and let TTk^i be the /-th coordinate of p^, 0<i^k<n — 1 and 
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1 ^ i, ^ ^ f^- To decide which one of the two corresponding perturbed coordinates is smaller, we 
define a predicate Smaller as follows: 

Smaller [tt i ; TTk^i) = true iff T^i,j{£) < T^k,i{£)- 

Due to SoS, we can neglect degeneracies, i.e., we have vrij(e) ^ 7r/;^;(e), and for this reason the pred- 
icate Smaller [jTi^j; T^k,i) = false if and only if vrij(e) > TTk,i{£)- The implementation of this predicate 
is fairly straightforward since we can compare the e-terms, e{i,i) and e{k,l), by comparing the 
defining index pairs (see Section 3.2, Lemma 3.2). 

Predicate 1 {Smaller) Assume the e-expansion e(z,j) defined as in Section 3.2 (2). With this, for 
indices 0<z,A;<n — f and 1 < j ,1 < d which satisfy (z, j) ^ (A;, /), the predicate Smaller [iTi^f^ '^k,i) 
can be implemented as follows. 

function Smaller {i^i^j]i^k,l) returns Boolean 
begin 

if TT^-^j 7^ Tik,l then 

return (vr^-^j < Tik,l) 
else \i i ^ k then 

return {i > k) 
else 

return (j < /) 

end 

Notice that, in this case, the coordinates tTj-j and TTk,i as well as their index pairs and (A;,/) 

have to be passed as arguments whenever predicate Smaller is called. This means that in popular 
programming languages, such as Pascal, the function heading would be something like 

FUNCTION smaller (i, j, k, 1, Pij , Pkl) : Boolean; 

but implementation details like this will be ignored in the remainder of the paper. Furthermore, 
notice that we have 

Smaller (iTi TTk i) = true iff det ( | ) < 0. 

\ T^k,i{£) J- / 

In Section 3.1, we will express more complicated predicates than just comparisons of coordinates 
by similar determinants. For matrices not exceeding a given size it is not difficult to specify the 
e-expansion e(z,j) such that all requirements discussed in Section 2 are satisfied. This will be 
done in Section 3.2. Finally, Section 3.3 extends the results to homogeneous coordinates. The 
procedures that implement the predicates will be developed in Section 4. 

3.1 Predicates Expressed by Determinants 

This section introduces the notion of orientation of a sequence of J + I points in E'^. With this 
concept, we will be able to give an implementation of the primitive operation for d -\- 1 points 
mentioned above. 
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The orientation of a sequence of points (pi^^pi^, . . . ,pij) in E'^ is either negative or positive — 
unless the d -\- 1 points lie in a common hyperplane, in which case the orientation is undefined. 
The exceptional case is a degeneracy that can be ignored if the points are perturbed. We define 
the orientation of a sequence recursively. It will be important that the orientation of a sequence 
depends only on the relative position of the points to each other and not on their absolute positions. 

If the dimension J = 1, then the orientation of {pio^Pi^) is positive if pi^ > pi^ and it is negative if 
Pio < Pii (compare with Figure 3-1, (a) and (b)). If J = 2, then [pi^^pi^ 7^82) has positive orientation 
if the three points define a left-turn in the plane, that is, pi^ lies to the left of the directed line that 
passes through pi^ and pi^ in this order. If [pi^ ^ pi^ ^ pi^) defines a right-turn, then its orientation 
is negative. Note that the orientation of [pi^^pi^^pi^) is the same as the orientation of {pi^,Pi2) 
as "seen from" pi^ . Indeed, the line through pi^ and pi^ can be identified with as soon as we 
choose a direction of the line. This direction is provided by the location of pi^ : It goes from left 
to right as seen from pi^ (see Figure 3-1, (c) and (d)). 



Pii 




Pn P^o P^o Pn P^o' 

(a) positive (b) negative (c) positive (d) negative 

Figure 3-1: The orientation of J -|- 1 points in dimension J, for d = 1,2. 

If J > 2, then the orientation of [pi^^pi^^ . . . ,Pi^) is the same as the orientation of {pi^, ■ ■ ■ ,Pi^) 
as seen from pi^. For example, (pi^^pi^^pi^^pi^) in E^ has positive orientation if pi^ observes 
{PinPi2 7Pi3) making a left-turn. In most situations where the concept of orientation is used, the 
interest is in the position of one point, pi^^ relative to d other points, pi^^pi^^ . . . ,pi^. We thus say 
that pi^ lies on the positive side of[pi^^ . . . ,Pi^) if (pi^^pi^^ . . . ,Pi^) has positive orientation, and pi^ 
lies on the negative side of [pi^^ . . . ,pi^) if [pi^^pi^^ . . . ,pi^) has negative orientation. 

To decide upon the orientation of a sequence of J -|- 1 points in E"^ ^ we use the matrix 



A 



(3-a) 



Lemma 3.1 The orientation of [pi^^pi^^ . . . ,Pi^) is positive if and only if sign(det A) = -|-1 and it 
is negative if and only if sign(det A) = —1. 

Notice that det A vanishes if and only if the d -\- 1 points are degenerate, that is, they lie in a 
common hyperplane — a case that can be neglected within the perturbed point set -P(e). Recall 
from linear algebra that the determinant of a matrix is multiplied by —1 if we exchange two rows. 
Thus, the orientation of a permutation of [pi^^pi^^ . . . ,Pi^) is the same as the orientation of the 
sequence itself if the number of transpositions is even, otherwise, its orientation is the opposite of 
the orientation of [pi^ , , . . . , Pi^)- 



Simulation of Simplicity 



8 



There are plenty of algorithms for point set problems which are based on computing the orientation 
of a sequence of points. Prime examples are the construction of convex hulls (see [PH77], [PS85], 
[Se81], [Se86], or [Ed87]), computing A-matrices as discussed in [GP83] and [Ed87], and finding 
convex subsets (see [CK80], [EG89], and [Ed87]). The remainder of this section considers the 
primitive operations required by the three-dimensional convex hull algorithm of Preparata and 
Hong which is described in [PH77], [PS85], and [Ed87]. 

The first step of the algorithm sorts the points in Xi-direction. To perform this step, it needs to 
compare the Xi-coordinates of two points, which can be done by computing the orientation of their 
orthogonal projections onto the Xi-axis. Second, it constructs the two-dimensional convex hull of 
the points projected onto the Xi^s-plane. Here, the primitive operation is to decide whether three 
points (in the Xi^s-plane) define a left-turn or a right-turn. Third, the algorithm constructs the 
three-dimensional convex hull by repeating the following operation: 

Given a plane pivoting about two extreme points pi^ and pi^ , find the point hit first 
by this plane. 

This operation can be reduced to a number of comparisons of the form: Given two points pi^ and 
Pig , which one is hit earlier by the pivoting plane? To perform such a comparison is equivalent 
to deciding on which side of the plane through p^g, pi^^ and pi^ point pi^ lies. This is the same 
as computing the orientation of (pi^^pi^^pi^^pi^). Thus, we see that the convex hull algorithm of 
Preparata and Hong requires three primitive operations all of which determine the orientation of 
point sequences. 



3.2 Choosing the Form of the Perturbation 

As explained in Section 3.1 the primitive operation that determines the orientation of a sequence 
of d-\-l points in d dimensions computes the sign of a determinant of a [d-\- l)-hy-[d -\- 1) matrix. 
SoS replaces the coordinates tTj-j in this matrix by entries of the form tTj-j -|-e(z, j). The determinant 
itself is then the sum of a finite number of terms, where each term is the product of d items and an 
item is either an original coordinate or an e{i^j). Thus, each term consists of a coefficient, which 
is the product of original coordinates, and a so-called e-product, a product of factors of the form 
e{i^j). The number of factors e{i,i) can be zero in which case the e-product is defined to be equal 
to 1. As mentioned in Section 2, it is irrelevant what exactly the definition of the e-expansion 
is as long as it satisfies certain requirements. The computational simulation is uneffected if we 
change the definition of the e-expansion within allowed limits. Even so, it is important to show 
that there is at least one e-expansion that satisfies the requirements. The existence of such an 
expansion implies the physical existence of an appropriately perturbed point set which is the only 
guarantee for the consistency of our method we have. 

We define 

for 0<z<n — l,l<j<J, and 6 > d, and show that this choice satisfies all the requirements 
of SoS. Notice that the amount of perturbation experienced by coordinate tTj-j is larger than the 
perturbation of TTk,i if and only if (z, j) -< (A;, /), that is, i < k ov i = k and j > I. Eurthermore, we 
have 

n ^(^^O = n e'''-'>e''''-'=e{Kl) (3-c) 

(t,3)<(k,l) (t,3)<(k,l) 
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if < e < 1. This is equivalent to stating that 2^'^~\ the exponent of e{k,l), is larger than the 
sum of the exponents of all e{i,i) with (z, j) -< (A;, /). It follows that it is sufficient to consider the 
sets of index pairs when we compare two e-products. Let ei and 62 be two different e-products 
and let X(ei) and 1{e2) be the two associated sets of index pairs. We call X(ei) smaller than 1{e2) 
if the set X(ei) — 1{e2) is empty or if (z, j) -< (A;, /), for (z, j) the largest index pair in X(ei) — 1{e2) 
and (A;, /) the largest index pair in 1{e2) — 1{ei). 

Lemma 3.2 Let Ci and C2 be two positive constants and let ei and 62 be two different e-products. 
Then Ci • ei > C2 • 62 for a small enough e if X(ei) is smaller than X(e2). 

Lemma 3.2 is an immediate consequence of (3-c) and the fact that a small enough e can compensate 
the influence of the constants Ci and C2. Notice that it is actually irrelevant which index pairs X(ei) 
and 1{e2) contain. The only thing of importance is the relative position of X(ei) and 1{e2) in the 
ordering of all sets of index pairs, where large index pairs are more significant in the comparison 
of sets than small index pairs. Observe also that Lemma 3.2 holds if we increase the value of 8 in 
the definition of the e-expansion. It turns out that this lemma is the crucial property that allows 
us to prove that -P(e), the perturbed point set, is simple and that the orientation of J + I points 
in P{e) can be computed efficiently. 



Lemma 3.3 The set P{e) is nondegenerate if e > is sufficiently small. 



Proof. To prove the assertion, we show that for no choice of J + I mutually distinct indices 
Zo, Zi, . . . , id, the determinant of the matrix 



/ 



2'o^ 



2'o^ 



vr,:, 1 + e tt,^^2 + e 



2'd ' 



2'd' 



2'o- 



2'1' 



2'd- 



1 \ 



1 / 



(3-d) 



is equal to zero. To see this, we assume w.l.o.g. that 0<io<ii<...<id<n — 1 and sort the 
terms of det A(e) in order of increasing exponents of e. Specifically, 



det A 



is the first term, and 

^_-^yd/2] _ ^2'1 *-'*-|-2'2 «-(d-l) + ... + 2'd «-l 

the last one. Each term is of the form for some constants b and c. Because we can assume that 
e > is arbitrarily small, the absolute value of the first term with nonzero coefficient b is bigger 
than the sum of all other terms. Furthermore, such a term always exists since (3-c) guarantees 
that no two terms of the determinant have the same exponent of e and thus such a term cannot 
cancel. For example, the coefficient of the last term is ( — l)r<i/2l Q and cannot be canceled by 
any other term. Consequently, det A(e) does not vanish. □ 

As pointed out in the proof of Lemma 3.3, the most significant term of the polynomial det A(e) 
is the determinant det A of the original coordinates. If the orientation of the original sequence 
{Pio 7 Pin ■ ■ ■ 7 Pid) defined, then this term is nonzero which implies that the orientation of the 
perturbed sequence is the same. This is reassuring since it shows that the perturbation does not 
change nondegenerate relations of the original point set. 
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The curious reader might wonder why the perturbation is defined in the pecuhar form given by 
the e-expansion (3-b). As mentioned before, there are many other choices that could be used, for 
example, 

e(z,j) = e 

is such a possibility. This e-expansion would also work but its implementation is slightly more 
difficult than that of (3-b) (compare with Section 4.2). On the other hand, many less "exotic" 
choices do not work. The remainder of this section illustrates this by considering two choices of 
e(z,j) which appear simpler than (3-b). The two choices are 

e{i^ j) = e''^~^^ and e{i^ j) = [i ■ 6 + j) ■ e. 

In both cases. Lemma 3.3 does not hold. The reason for the failure is that both expansions do 
not satisfy (3-c) and thus possibly lead to cancellations of e-terms in det A(e). Such cancellations 
occur for example if all J -|- 1 points of the sequence coincide with the origin. In this case the 
matrix A(e) equals 

/ e(zo, 1) e(zo,2) • • • e{io,d) 1 \ 
e(zi,l) e(zi,2) ••• e{ii,d) 1 

V e{id, 1) e{id,2) ■ ■ ■ e{id,d) 1 / 

If we define e{i,j) = e^'^~^^ , then the second column is equal to e times the first column which 
implies that det A(e) = if J > 2. If e(z, j ) = {i • 6 -\- j) ■ then the sum of the first and the third 
columns equals twice the second column; hence, det A(e) = if J > 3. 



3.3 Homogeneous Coordinates 

When we develop the primitive procedures for computing the orientation of d-\-l points in Section 4, 
we represent a point by its homogeneous coordinates. This representation is slightly more general 
than ordinary Cartesian coordinates (it can also represent points at infinity) and leads to a slightly 
more uniform procedural treatment. 

Let p be a point in E'^ and let (vr^, tt^, . . . , tt^) be its sequence of Cartesian coordinates. Point p 
has d -\- 1 homogeneous coordinates 

(TTi ,7r2 ,...,7r^ ,7r^+ij 

such that 

H 

vrf = for l<i<d. 

■ 

Thus, p is I/tt^j^-^ times the point whose Cartesian coordinates are equal to the first d homogeneous 
coordinates of p. Notice that the homogeneous coordinates of p are not unique; we still represent 
the same point p if we multiply each coordinate by the same nonzero scalar. If we decrease the 
absolute value of vr^^^ without changing the other homogeneous coordinates, then p moves away 
from the origin on a straight line and it reaches "infinity" when vr^^^ becomes 0. Indeed, p is "at 
infinity" if and only if tt^_^i = 0. Using homogeneous coordinates, it is not allowed to have all 1 
coordinates are equal to — in this event p is not defined. 
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We next extend Lemma 3.1 to homogeneous coordinates, that is, we characterize the orientation 
of a sequence of J + 1 points [pi^^pi^ , . . . 

in terms of their homogeneous coordinates. The orientation of a sequence of J + 1 points is not 
defined if any of the points Res at infinity. In fact, it is not possible to generahze the notion of 
orientation to points at infinity without changing our interpretation of a point at infinity. For 
example, consider a sequence S of d finite points and one point p = (tt^ , tt^, . . . , ir^; 0) at infinity. 
We can think of p as the hmit of points 

= (Trf ,vrf , . . . ,7rf ; e), 

when e > goes to zero, but as well, we can think of p as the limit of these points if e is negative 
and approaches zero. If we replace p by p(e) with e small enough, then e > and e < lead to 
different orientations. We thus restrict our discussion of orientation to finite points. Define 



A = 


/ 


<1 




■ it" 




(3 




v 


<1 






) 





If irf^ = 1, for < V < d, then A is the same as the matrix A used in Lemma 3.1. Otherwise, 
we can multiply the rows such that ir^ = 1. The sign of det A changes if we multiply a row 
with a negative number, which implies the following result: 

Lemma 3.4 Let {p^,,p^,,. . . be a sequence of points with p,^ = [tt^^^-^, irf^^^ . . . , irf^^^, T^t,d+i) 
and TT^ ^ 0. Their orientation is positive if sign(det A) = Ylt=o^'^S^{'^il d+i)^ negative if 
sign(det A) = — Ylt=o ^'^S^i'^il d+i)^ undefined if det A = 0. 

In contrast to Cartesian coordinates, a point is now represented by J + 1 coordinates which makes 
it necessary to choose 6 > d -\- 1 when defining the e-expansion e{i,i) in (3-b). With this, it is 
easy to prove that determinants cannot vanish which implies that Lemma 3.3 holds also for the 
new setting using homogeneous coordinates. 

4 Implementing a Predicate 

This section presents the actual implementation of a geometric predicate using SoS. The chosen 
predicate determines the orientation of a sequence of points, as defined in Section 3. Its imple- 
mentation will be based on the e-expansion specified in Section 3.2 (3-b) and on the fact that the 
orientation can be found by evaluating the sign of a determinant as stated in Sections 3.1 and 3.3. 
The crux of the implementation is that this determinant is a polynomial in e. The computation 
of the sign of such a polynomial is discussed in Section 4.1. The coefficients of the polynomial 
turn out to be subdeterminants of the original matrix. Based on this observation. Section 4.2 
gives an algorithm that generates these subdeterminants in sequence of decreasing significance by 
employing a special encoding scheme. Finally, in Section 4.3 we will briefiy address the problem 
of sign computation of integer determinants in general. 
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In Sections 3.1 and 3.3 we defined tlie "orientation" of a sequence of points in J-dimensional 
Euclidean space given by Cartesian and homogeneous coordinates. We now formally develop the 
corresponding predicate that uses perturbation in the sense of SoS. In the Cartesian case each 
point is given by its d coordinates 

whereas in the case of homogeneous coordinates a point is represented by a ( J + l)-tupel 

Piy = (tTj.,!, . . . , T^,y,d; T^,y,d+l)- 

Let 

P = {po, . . .,Pn-l} 

be a set of n points in E'^, and denote by 

^(e) = {Po(e),...,Pn-i(e)} 

its perturbed version using the e-expansion of Section 3.2 (3-b), assuming 6 large enough such 
that Lemma 3.2 is valid. Now define for d -\- 1 points with distinct indices Zq, Zi, . . . , Zd, all in the 
range from through n — 1, 

Positive d{pio^ ■ ■ ■ ^Pia) = true iff the orientation of (^^^(e), . . . ,pi^(e)) is positive. 

Degenerate cases can be neglected because we simulate simplicity. From Lemma 3.1 it follows that 
Positive d is equivalent to the test whether or not 

sign(det A(e)) = +1, 

with A(e) denoting the the corresponding matrix of the perturbed Cartesian coordinates as in 
(3-d). In the homogeneous case (see Lemma 3.4) we have to check whether or not 

d 

sign(det A(e)) = sign(7r,^,d+i(e)). 

i/=0 

Here, A(e) denotes the perturbed version of matrix A in (3-e), whose rows are formed by the 
homogeneous coordinates of the points involved; that is, 

/ TT^o,! + e(zo, 1) 7rjo,2 + e(«o,2) ••• ir.^^d+i + £{io,d + 1) \ 
TTjj,! + e(zi, 1) 7rjj,2 + e(«i,2) ••• vr^j^d+i + e(zi, J + 1) 
^y^)= : : •. : 

\ + £{id,l) 7rj^,2 + e(«d,2) ••• vr^^^d+i + e(zd, J + 1) / 



At first sight, the development of such an e-determinant seems to be a painful exercise. Yet, it will 
turn out that it is not that hard and can be achieved in an algorithmically clean way. Anyway, to 
begin with something easy, consider 

Hpt A (^) - Hpt ( ""''1 + ""''2 + 2) ^ 
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Let e((zi,ji), . . . , [ik^jk)) = Ylt=i ^(^i'jii')? call it a k- f old e- product; e() = 1 is called the 0-fold 
e-product. Furthermore assume i < j. When we now develop the determinant we get 

detA,{e)= +f^''^ ""^'O-^O- 

- TTj^i ■ e{i, 2) + 7rj,2 • e(«, 1) + (4-a) 
+7r,,i-e(j, 2) + l-e((j, 2),(^,l))- 
-7r,■2•e0• l)-l-e((j, l),(z,2)), 

where the terms are already sorted by increasing powers of e. Note again that the first coefficient 
corresponds to the "unperturbed" determinant, i.e., A2, whose evaluation would be part of any 
implementation of the predicate — of course, followed by the more or less awkward handling of 
all possible degeneracies. Observe also that the coefficient of the fifth term is a constant, namely 
+ 1. Thus, the last two terms have no infiuence on the sign of det A2(e). Therefore, the number 
of relevant terms of the e-polynomial det A2(e) is only 5, rather than 7 which is the total number 
of terms. 

It is convenient to assume io < . . . < id (compare with (4-a)). This assumption together with 
Lemma 3.2 implies that the sign of det A(e) and det A(e) can be computed without any further 
knowledge of the values of the indices. Clearly, this is not the case in general but can always 
be achieved by appropriate row exchanges in A(e) or A(e) — recall that each exchange changes 
the sign of the determinant. For this, assume a procedure Sortd+i{{io^ . . . , Zd), (zq, . . . , z'^), s') that 
returns for a given sequence of J -|- 1 indices (zq, . . . ,id) the sorted sequence (zq, . . . , z'^). Addition- 
ally, Sortd+i returns s' which is set to the number of exchanges used. We can now implement 
predicate Positive d using two operations, SignDetA and SignDetA^ that compute the sign of the 
e-polynomials det A(e) and det A(e) assuming io < . . . < id- Both functions will be discussed in 
Section 4.1. 



Predicate 2 (Positive) Let pi^, . . . ,pij be 1 points in E'^ given in Cartesian or homogeneous 
coordinates with distinct indices all between and n — 1. Then the following pseudocode is an 
implementation of the predicate Positivcd- 
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function Positive^ {pi^, . . . ,pi^) returns Boolean 

local ig, . . . , i^, d' , s' , v 

begin 

SortdJ,i{{iQ, id), i'^), s') 

if Cartesian coordinates then 

d' ^ SignDetAd+i 



else 



d' ^ SignDetAd+i 



if odd(s') then d' i d' 

if Cartesian coordinates then 

return {d' = +1) 
else 

return {d' = flLo sign(7r,^,d+i(e))) 

end 



The problem is now to give efficient implementations lor the two lunctions SignDetAd+i and 
SignDetAd+i- We leel that it is important to stress that "efficiency" is meant in a practical sense 
— in theory it can be done in constant time anyway, assuming J is a constant. 



4.1 The Sign of a Perturbed Determinant 

We now illustrate the implementation of SoS on the bottommost programming level by imple- 
menting function SignDetAo, which returns the sign of a D-hy-D e-determinant det AoiA for 
any given D; primitive SignDetAo can be treated in the same way. To appreciate the significance 
of a (practically) efficient implementation of SignDetAo we point out that this is in fact the major 
part of SoS, at least when applied to the predicate described above. Provided that Zq < • • • < zd, 
we will show that it is possible without great effort to generate the sequence of the coefficients of 
det AoiA in decreasing order of significance. Since e can be assumed to be sufficiently small (but 
positive), the sign of the e-polynomial is therefore equivalent to the sign of the first nonvanishing 
coefficient. 

Using simple rules for evaluating a determinant as exemplified for det A2(e) in (4-a), the coefficient 
of every term in det AoiA is a subdeterminant of the "unperturbed" matrix A^- Here, a single 
entry is called a 1-by-l subdeterminant and, by definition, the 0-by-O subdeterminant is equal to 
1. To tell the whole truth, we must mention that each coefficient in effect is a subdeterminant 
together with a certain sign, that is, multiplied by either +1 or —1. We will see in Section 4.2 
how to decide whether +1 or —1 applies. To continue our discussion, we need a few notations. 
We say that the (t + l)-st coefficient in order of decreasing significance, denoted by det Af^"^^, is 
the cof actor of depth t of matrix AoiA- Note that this coefficient already includes its proper sign. 
Thus, det = -j-detAu- The size of the corresponding matrix (i.e., the number of rows or 

columns) is denoted by kt = k{M^°). These definitions are illustrated in Table 1 which shows all 
significant terms of det A2(e). In the column with heading St we display the e-product associated 
with the cofactor of depth t. The column Vt will be explained later. 
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t 


kt ■ kt 




det 2 




u 


2 • 2 


[3,3;3J 


+ det 


e() 






V VTj- 1 TTj- 2 y 


1 


1 • 1 


[2, 3; 3] 


— det (vTj^i) = — VTj^i 


e(^,2) 


2 


1 • 1 


[1,3; 3] 


+ det (7rj^2) = +7rj,2 


e(^l) 


3 


1 • 1 


[2, 2; 3] 


+ det (vr^^i) = +7r8,i 


e(j,2) 


4 


0-0 


[1,2; 3] 


+ det() = +1 


e((J,2),(^,l)) 



Table 4-i: The 5 relevant terms of det A2(e). 



This leads to the pseudocode implementation of SignDetAo shown below. It assumes that Zq < 
. . . < ijj and that the sequence of subdeterminants, sorted by increasing depth, is known. The 
code also requires a function SignDetj^{^) that calculates the sign of det $ for a k-hy-k matrix $. 
The authors have not been able to find an alternative way to determine the sign but to compute 
the actual determinant. Unfortunately, computing the (exact) determinant of a matrix of integers 
demands the use of long integer arithmetic. More about that in Section 4.3. 



function SignDetAu (A/j) returns +1 or —1 

local (7, kt, t 

begin 

t ^ -1 
repeat 

t '-t+l 

kt ^ k{M^° ) 

until (7 7^ 
return a 
end 



scans 



Function SignDetAu 

pseudocode, "fct ^ k{M^°y and "ci SignDetf.^(M^ 
could be implemented CASE-statement. For D = 
shown below. 



through the table of relevant subdeterminants. Two lines of the 



^)", indicate table lookups. In Pascal this 
2, it would consist of 5 different 



CASE t OF 

1 
2 
3 
4 

EID; 



SignDet2 (Pil ,Pi2 ,Pj 1 ,Pj2) ; 
-Sign (Pjl) 
Sign (Pj2) 
Sign (Pil) 

1; 



If the depth counter is of no interest, one can even unwind the loop and come up with the following 
code. 
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FUICTIOI SignDetDelta2 (Pil, Pi2, Pjl, Pj2) : Integer; 
BEGII 

SignDetDelta2 := SignDet2 (Pil, Pi2, Pjl, Pj2); 
IF SignDetDelta2 <> THEI goto 999; 
SignDetDelta2 := -Sign (Pjl); 
IF SignDetDelta2 <> THEI goto 999; 

SignDetDelta2 := 1; 
999: (* exit *) 



To give more insight into the computation of the terms of det Aj:)(e) in the order of decreasing 
significance, we now consider the three-dimensional case, that is. 



This polynomial has a total of 34 terms. However, only 15 of them are relevant, and those are 
listened in Table 2. There are two reasons why we only need to test 15 coefficients out of a total 
of 34. One is that the coefficient of e((A;,3), (i,2), (z, 1)) is equal to +1 which is nonzero; we can 
therefore stop there and consider no further terms. The other reason is that certain coefficients 
occur more than once, that is, with different e-products. For example. 



Clearly, there is no need to test — vr^^s ^ 0, since at this depth, -\-Trk,3 = is already known; 
otherwise, the sign determination would have stopped immediately after testing the coefficient of 



4.2 Generating the Sequence of Significant Coefficients 

The properly sorted sequences of e-terms of the polynomials det A2(e) and det A3(e) are apparently 
very regular. In the following, this regularity will be worked out and exploited by an algorithm 
that automatically generates the correct sequence of e-terms. This procedure can be embedded 
in an implementation of the function SignDetAo that computes the sign of det Aj:)(e). We agree 
that a procedure that generates each term of det Aj:)(e) by collecting the proper rows and columns 
of the original matrix is, in a practical sense, much slower than a straight-line program that scans 
through a fixed sequence of submatrices. However, in higher dimensions the former might be the 
better strategy, since the likelihood of det = for all r with < r < t decreases very fast 

as t increases. Let alone the fact that the tables of relevant terms for det Aj:)(e) becomes rather 
long for large D. The algorithm to be described can also be used for automatic generation of such 
tables and even for the automatic generation of codes implementing them. 

We now discuss in detail how we can extract the individual terms of the polynomial det Aj:)(e). 
Recall that a term is of the form 6 • e^, where b is called the coefficient and is the e-product of 
the term. If = e((zi, ji), . . . , (ik^jk)) (so it is a A;-fold e-product), then we call e{i^^j^) active^ 
for 1 < L < k. Given the e-product of a term we can extract the coefficient b from the given 



EID; 




det A3(e) = • • • + 7rfc,3 • e((j, 2), (z, 1)) 7rk,3 ■ e((j, 1), («, 2)) • • • 



(4-b) 



£((j,2),(z,l)). 
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t 






det Mf 3 







3 • 3 


[4, 4, 4; 4] 


/ T^t,! 1^1,2 1^1,3 \ 
+ det TTj^i 7rj^2 TTj^s 

\ TTfc^i 7rfc^2 TTfc^s y 




1 


2 • 2 


[3, 4, 4; 4] 


+ det f ''^■'^ ''^■'^ ^ 


e(i,3) 


2 


2 • 2 


[2, 4, 4; 4] 


- det f ""^'^ ""^'^ ] 


e(i,2) 


3 


2 • 2 


\} 4 4- 41 

[-■-7 ^7 ^5 


+ det ( ""^'^ ""^'^ ] 




4 


2 • 2 


[3 3 4- 41 

[O, O, 4:, 4:J 


- det f '''■'^ ""''^ ] 


f( i 3^1 


5 


1 • 1 


[2, 3, 4; 4] 


+ det (vrfc^i) = +TTk,i 


e((J,3),(^,2)) 


6 


1 • 1 


[1,3, 4; 4] 


- det (7rfc,2) = -vrfc,2 


e((J,3),(^,l)) 


7 


2 • 2 


[2, 2, 4; 4] 


+ det f '''■'^ ""^'^ ] 

\ T^k,l T^k,3 J 


e(j,2) 


8 


1 • 1 


[1,2, 4; 4] 


+ det (vTfc^s) = +TTk,3 


e((J,2),(^,l)) 


Q 


2 • 2 


[1 1 4-41 

[±, ±, 4:, 4:J 


- det f '''■'2 ""''^ ] 






2 • 2 


To O O . /1 1 

[3,3,3;4J 


+ det 

V VTj- 1 TTj- 2 y 


S[K, 6) 


11 


1 • 1 


[2, 3, 3; 4] 


— det (vTj^i) = — VTj^i 


e((A;,3),(^,2)) 


12 


1 • 1 


[1,3, 3; 4] 


+ det (7rj^2) = +7rj,2 


e((A;,3),(^,l)) 


13 


1 • 1 


[2, 2, 3; 4] 


+ det (vr^^i) = +TT,^i 


e((A;,3),(j,2)) 


14 


• 


[1,2, 3; 4] 


+ det() = +1 


e((A;,3),(J,2),(^,l)) 



Table 4-ii: The 15 relevant terms of det A3(e). 

matrix by crossing out all rows and columns that contain an active e{i^^j^). In order to avoid 
extensive double indexing and index inversions, we assume that the points whose coordinates are 
the entries in the D rows of the matrix have indices 1 through D. This allows us to ignore 
the difference between a point index and the corresponding row index. Indeed, this assumption 
is no loss of generality since the only property used in computing the sign of det Aj:)(e) is that 
the point indices are sorted and therefore the actual values are irrelevant. With this assumption, 
e(zt,jt) is in the z^-th row and the jVth column and we cross out rows Zi,Z2, and columns 

ji, j2, • • • This leaves a [D — k)-hy-[D — k) submatrix. Table 2 illustrates these definitions for 
D = 3. If 6 • is the term of depth t, then the notation in Table 2 is such that h = det Mf^° , 
= Ct, and kt is the number of rows (or columns) of . 

Note that we did not yet specify how we can decide whether 6 is —1 or +1 times the determinant 
of the submatrix. We now describe a rule that is based on the number of transpositions needed to 
sort a certain permutation. For row i, 1 < i < -D, let be the column such that e{L^j^) is active 
in the term that we currently consider. By definition of a determinant there can be at most one 
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such column but it could very well be that there is no such column. In this case we choose such 
that TTtj-^ belongs to the main diagonal of the submatrix that was obtained after crossing out rows 
and columns as described above. If the number of exchanges needed to sort (ji, J2, ■ ■ ■ ^jo) is odd, 
then b = det Af^"^^ is —1 times the determinant of the submatrix, otherwise, it is +1 times this 
determinant. 

Interestingly, the number of exchanges needed to sort the sequence (ji, J2, ■ ■ ■ ^jo) is even if and 
only if ii_ -\- ji^ is odd for an even number of pairs {it^ji,), 1 < l < k. To see this notice that the 
total number of pairs (k, Jk) with k, + odd is even since 

D D 

Now observe that (ji, J2, • • • ^jo) can be sorted using only exchanges of adjacent columns, that is, 
of integers that differ by one. Note also that we can dispense with all exchanges between two 
columns where both contain an active e(z, j) or both do not. Thus, every exchange of two columns 
increases or decreases the number of pairs with odd by one, which implies the claim. 

This property will be used in the algorithm that computes the proper sign. 

The key observation that allows us to automatically generate the relevant terms of det Aj:)(e) 
is that e((zi, ji), . . . , (ik^jk)) is the e-product of a relevant term if and only if ii < . . . < ik and 
ji < . . . < jk- In other words, the e{i^^j^) go monotonically from the left top to the right bottom 
of the matrix. To see this take an e-product that does not satisfy this condition and consider 
the e-product defined by the same 2k indices that is obtained by matching the smallest with 
the smallest j^, the two second smallest indices, etc. This new e-product is more significant than 
the old one since the exponent of e it defines is smaller than the exponent of the old e-product. 
Furthermore, the coefficients that correspond to the two e-products have the same absolute value, 
namely the determinant of the submatrix obtained by crossing out rows and columns j^, for 
1< L<k. 

The algorithm that generates the e-products and their corresponding coefficients uses a vector 

V = [vi,. . . ,vd;vd+i] 

where each Vi is an integer between 1 and D -\- 1 and Vi corresponds to the z-th row of det Aj:)(e); 
V£)+i is set equal to D -|- 1 and is used only for convenience. The interpretation of v is as follows: 
To encode the e-product e((zi, ji), . . . , (ik^jk)) we set Vi^ = for 1 < t < k. For every z such that 
the z-th row does not contain an active e{i^j)^ we define Vi = Vi^ with z^ the smallest integer in 
{zi, . . . , Zfc, -D -|- 1} that is larger than z. Thus, in v implies that e(/€, v^) is active if and only if 
< For example v = [3, 4, 4; 4] implies that the e-product of the encoded term is e(l,3). 

Other examples can be found in Table 2, which gives the vectors of all relevant terms in det A3(e). 

The next problem that we face is how to generate the terms of det A/^ in the correct order, that 
is, in the order of decreasing significance. Here we use the fact that v = [ui, . . . , u^; t^D-i-i] encodes 
a more significant term than v' = [uj, . . . , u^; t^/^+i] if and only if > v'^ for j the largest index, 
such that ^ v'y This implies that v = [D -\- 1, . . . , D -\- 1] D -\- 1] encodes the most significant 
term and, indeed, it encodes e() = 1, whose coefficient is the determinant of the entire original 
matrix. It is now easy to write a function that computes for a given vector its successor. 
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function Next-V (v) returns Vector 

local i, K 

begin 

L ^ 1 

while f t = 1 do i ^ i + 1 

V, ^ v,-l 

for K <— i — 1 down to 1 do ^ 
return v 
end 

The alert reader will have noticed that the above function returns an "illegal" vector if the input 
vector is [1, . . . , 1; D + 1] which is not a problem because the determinant evaluation is such that 
already [1, 2, . . . , D; D + 1] encodes a nonzero coefficient thus there is no reason to call Next-V 
again. 

After initializing u to [D + 1, . . . , D + 1; D + 1], successive calls to Next-V give the desired sequence 
of vectors. It remains to be shown how the coefficient of the encoded term can be computed. The 
procedure below decodes v and returns the submatrix M obtained after deleting the proper rows 
and columns from A/^. It also returns s equal to —1 or +1 depending on whether the coefficient 
equals — det Af or + det Af , and returns k which is equal to the number of rows (or columns for 
that matter) of M . 

procedure Matrix (v, s, k, M) 
global Ac, D 
local i 
begin 

M ^ Ad 
k ^ D 
+1 

for i ^ 1 to d do 
if ft < then 

{in this case e(/,,ft) is active} 
if odd(/, + ft) then s <— —s 
delete row i from M 
delete column from M 
k ^ k-1 

end 

We can now modify the code of SignDetAo by replacing the table lookup by appropriate calls to 
Next-V and Matrix. With additional modifications the same algorithm can be used to generate the 
table of relevant terms in det Aj:)(e) or even to generate the corresponding code for SignDetAo for 
any D. Note that, in the latter case, the loop in SignDetAo is to be repeated only until kf = 0, 
since in "generating mode" the values of the determinants are not computed and thus there is no 
natural abortion of the cycle of calls. The result for D = 4 can be seen in Table 6 in the Appendix. 

A nice feature of the above algorithms is that we only need to change the initialization of v 
to [D, . . . , D; D] to get an implementation for SignDetAo which computes the sign of the e- 
polynomial detA(e). For this case, the loop over all relevant terms has to be repeated until 
either the corresponding cofactor is nonzero, or, if we are in "generating mode," until kf = 1. 
See Tables 3, 4 and 5 in the Appendix for the relevant terms of det Aoie) for D = 2,3,4. It 
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seems worthwhile to mention that Cartesian coordinates should be used whenever possible. This 
reduces the problem roughly by one "dimension" compared to the homogeneous case (compare 
e.g., Tables 2 and 4). 

The presented e-polynomials det Aj:)(e) and det Aj:)(e) illustrate that the computational overhead 
caused by SoS is acceptably small. One has to keep in mind that the most significant term of these 
e-determinants corresponds to the original determinant which expresses the primitive. So, there is 
no way around the evaluation of the sign of this determinant for any implementation. If the input 
data is nondegenerated the cost of SoS is obviously zero and, in general, it is rather unlikely that 
the polynomials have to be evaluated down to large depths. Indeed, the largest depth or the sum 
of all depths that occurs in a computation can be used as a measure for the degree of degeneracy 
of the input data. 

By evaluating the subdeterminants we systematically take care of all possible degenerate cases. 
Take for example the evaluation of det A3(e). Different cases can be distinguished by looking at 
the largest depth tmax reached during the computations. This tmax can be 0, I, 2, 3, or 4 and the 
corresponding degeneracy is as follows (compare with Table 4 in the Appendix): 

The three points p^, and pk are in general position. 

The three points are coUinear but pj ^ pk and the line containing the three 
points is not vertical. 

The three points lie on a common vertical line but p^ ^ pk- 
Point pj coincides with p^, but not with p^, and the line through pi and p^ is 
not vertical. 

All three points lie on a common vertical line and, p^ = pk- 

It would be interesting to see this somewhat unnatural case analysis in greater detail since it gives 
a nonobvious breakdown into degenerate cases that has curious properties. 

This discussion completes the implementation of SoS with respect to the predicate Positive d for 
point sets in E"^ . We considered both the Cartesian and the homogeneous case. The key was to 
find a method that generates the proper sequence of relevant terms of det Aj:)(e) and det Aj:)(e) 
ordered by decreasing significance. With this, the implementation of the functions SignDetAo and 
SignDetAo was easy. We will see in Section 5 that both functions can also be used to implement 
other predicates. 

4.3 Remarks on the Sign Computation of Determinants 

In the previous sections we reduced all computations to a sequence of sign evaluations of determi- 
nants. In the primitives discussed in this paper, the matrices are at most of size [d-\- 2)-hy-[d-\- 2), 
d the dimension of the space, and all elements are assumed to be integers. Theoretically, the sign 
of such a determinant can be determined in constant time if we assume that J is a constant. This 
assumption is indeed fair since SoS is intended primarily for low-dimensional geometric computa- 
tions. In practice, however, it is important to optimize the sign computation since it will be in 
the innermost loop of every program that uses SoS — which does not mean that this issue is less 
important for programs not employing SoS. We remark on a few methods that can be used to get 
speed in these computations. 

One important condition that we have to meet is that the sign of the determinant has to be 



^max • 

^max 1 • 

/ — 2- 

''max 

/ — 3- 

''max ^* 
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computed exactly — we cannot tolerate a +1 lor a 0, etc. Assuming that the coordinates or 
parameters are integers, we can either use long integer arithmetic or modular arithmetic based 
on the Chinese remainder theorem. For details on both methods reler to [Kn69]. II we actually 
compute the determinant in order to find its sign — and no method is known to the authors 
that avoids the actual computation of the determinant — we have to be prepared to deal with 
numbers of absolute size at least /i^, where // denotes maximum absolute value of any data item 
and D denotes the largest size of matrices we work with. To see this, just take the D-hy-D matrix 
whose entries are all zero except for the ones in the main diagonal, where they are equal to the 
determinant of this matrix is /i^. An upper bound on the absolute value of the determinants is 
given by a well known theorem of Hadamard that states that 



D 



|det(7ri..._D,i..._D)| < Yi 



\ 



Among other things this upper bound on the absolute value of a determinant gives us an upper 
bound on the number of computer words needed for the computation if we use long integer 
arithmetic. 

Without any hardware support long integer arithmetic is very time consuming, which might moti- 
vate us to resort to the use of approximation methods. Any computation of the determinant using 
fioating-point arithmetic of bounded length is such an approximation. Floating-point arithmetic 
is usually rather fast since it enjoys the needed hardware support on most of today's computers. 
If the value that we get is sufficiently far from zero, we can be sure that the correct value is dif- 
ferent from zero and lies on the same side of zero. But how can we quantify "sufficiently far from 
zero"? In any case, we could now use Gaussian elimination (see e.g., [GVL83]) which takes 0{D^) 
time or asymptotically faster methods based on matrix multiplication as described for instance 
in [AIIU74]. We do not believe that the latter methods could be of any practical use, though. 
However, if the value that we get is suspiciously close to zero, we have to use some other method 
to determine the sign of the determinant. 

Finally, we would like to mention that the determinant of a D-hy-D matrix can be expressed 
in terms of subdeterminants, and that some of these subdeterminants might later appear again 
when the evaluation of det Aj:)(e) or det Aj:)(e) proceeds. It is conceivable that the values of such 
subdeterminants are saved and used again when needed. Fven so, we do not believe that such a 
method could lead to significant savings since we expect that on the average only very few terms 
of the e-determinants are needed. 



5 Further Applications of SoS for Determinants 

In this section we demonstrate that the algorithmic solution to many geometric problems can be 
based on primitive operations that compute the sign of determinants. Those include problems 
that deal with objects different from points. There are two major reasons why determinants are 
useful beyond problems for points. One is that more complicated geometric objects are often 
given by a finite set or sequence of points. Fxamples are line segments given by two points and 
triangles specified by three points. This will be illustrated in Section 5.1, which revisits the Parity 
Algorithm discussed in the Introduction. The other reason (and this is the more profound although 
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less obvious of the two) is that other objects can be thought of as points in a different space. Take, 
for example, a hyperplane in d dimensions. It can be specified by a linear relation of the form 

f]ixi + f]2X2 H h r]dXd + Vd+i = 0. 

Multiplying the above relation with a nonzero constant does not change the hyperplane. This 
suggests that we think of the hyperplane as the point with homogeneous coordinates 

in d dimensions. This view of hyperplanes will be discussed in more detail in Sections 5.2 and 5.3. 
Of course, an n-gon specified by a sequence of n points in the plane can be interpreted as a point 
too — in this case it is a point in 2n dimensions. However, in contrast to the former case, this 
view is not likely to lead to any useful application of determinants since it becomes increasingly 
expensive to compute them as the size of the matrix increases. Finally, Section 5.4 shows that 
even nonlinear geometric objects such as circles and spheres can profitably be interpreted as points 
in low dimensions as well. 

By no means do we believe that the list of applications for primitives concerning the sign of 
determinants, as presented in this paper, is exhaustive. In fact, because of the versatility of 
determinants, an enumeration of their applications in geometric computation is far beyond the 
scope of this paper. We agree though that such an enumeration is a challenging task. 

5.1 Point-in-Polygon Test 

Recall the Parity Algorithm for the point-in-polygon problem sketched in the Introduction. In 
order to test whether a given point p lies inside a simple polygon P, the algorithm intersects the 
horizontal half-line r, whose left endpoint is p, with all edges of polygon P. If the number of edges 
intersecting r is odd, then p lies inside P, and if this number is even, p lies outside. The subtlety 
of this algorithm lies in the treatment of special cases since the above characterization holds, in 
general, only if we introduce certain artificial counting mechanisms whenever r contains a vertex 
or even an entire edge of P. In this section we show that the test whether or not an edge intersects 
the horizontal half-line r can be reduced to computing the signs of certain determinants. SoS is 
then used to simulate a perturbation of the point and the polygon which removes all degeneracies. 
The algorithm assumes that P is given by a sequences of vertices (ui, ^2, . . . , u„) and that all 
coordinates including those of p are integers. 

We consider now the problem to test whether r intersects an edge e of P given by its two endpoints. 
Let u = [vi^V2) and w = (c<;i,c<;2) be the two endpoints and recall that p = (7ri,7r2) is the left 
endpoint of r. Because of SoS we can assume that p are not coUinear and that no two of 

the three points lie on a common horizontal line. Note first that r and e intersect only if the 
second coordinate oi p lies between the second coordinates of u and w. Assume V2 < i02. If indeed 
V2 < TT2 < ^2 then r n e 7^ if and only if (u, w^p) defines a left-turn (see Figure 5-1). 

It is now not very difficult to develop this case analysis into a predicate that tests for intersection. 
To perturb the points we use the same e-expansion as described in Section 3.2, that is, we replace 
V, = {Pt,i,Pt,2) by ^^(e) = {i^t,i{£), i^t,2{^)) where Pt,j{e) = v,^j + e{ij) with e{ij) as in (3-b). For a 
uniform treatment we define p = Vq = (z^o,i, 1^0,2) and write the predicate for arbitrary three vertices 
rather than for Vq and two successive vertices of P. 
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w 




(a) (b) (c) 

Figure 5-1: The three cases to consider for r D e using SoS. 

In (a), r and e do not intersect since the second coordinate of p does not lie between 
those of u and w. In (c), they do not intersect since (u,w,p) is a right-turn. 



Predicate 3 (IntersectHalfLine) Let u^, Uj, and Ufc be three vertices with pairwise different indices 

< hj^k < n. The following pseudocode returns true if the edge from Vj[e) to Ufc(e) intersects 
the horizontal half-line whose left endpoint is given by t^i(e), and false otherwise. 

function IntersectHalfLine (vi;vj,Vk) returns Boolean 

local i' , j', k' , s' , d' 

begin 

W.l.o.g. assume SmaUer{vj^2]i^k,2)- 
if SmaUer{vj^2] i^i,2) A SmaUer{vi^2] ^k,2) then 
5ort3((^,J,'A;),(V,/,A;'),5') 

^ //j',i(e) Vt',2{£) 1 
vy,i{£) ^i',2{£) 1 
V Vk',i{£) Vk'ai^) 1 

if odd(s') then d' < d' 

return {d' = 
else 

return false 

end 

A few remarks are in order. When the above function is applied to the point-in-polygon problem, 
z = always holds. Thus, the sorting of A;) can be reduced to a single comparison between 

1 and k. Furthermore, to avoid all degeneracies for the point-in-polygon test, it is sufficient to 
perturb only the point p = Vq. Indeed, if 



d' <— SignDetAs 



det 



Z/j- 1 Z/j- 2 f 1=0, 

\ l^k,! '^k,2 



then we necessarily have Pj^2 = i^k,2 7^ ^i,2{£) and therefore, the determinant does not even get 
evaluated. The savings that one gets this way are only nominal which we interpret as an argument 
for the efficiency of our general method. 

The remainder of this section is used to comment on what happens if the test point p lies on the 
boundary of the polygon P. If we use the above primitive as is, SoS will neglect this special case 
and find that p lies on either side of P's boundary. The decision depends on the relative positions 
of p and the vertices of P, and we might as well assume that it is arbitrary although consistent. 
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Such a decision may or may not be desirable. If it is not acceptable, one could test whether or 
not p lies on the boundary of P before running the Parity Algorithm with SoS. Once more this 
test can be reduced to computing signs of determinants. 

5.2 Hyperplanes in Euclidean Space 

Algorithms for hyperplanes play a central role in computational geometry. This becomes obvi- 
ous when one thinks of the importance of problems such as linear programming, computing the 
intersection of half-spaces, and constructing arrangements of hyperplanes (see [PS85] and [Ed87] 
for further details and references). The goal of this section is to demonstrate how the techniques 
of Section 4 can be used to implement a typical primitive operation needed in those algorithms. 
This will open up an entire class of problems to the use of SoS. The main tool that lets us exploit 
the techniques of Section 4 when we handle hyperplanes is a duality transformation that maps 
hyperplanes to points and vice versa. In essence, this transform is nothing but a reinterpretation 
of what hyperplanes and points are. 

In this section we assume that a hyperplane h in E'^ is specified by its nonzero normal vector 
a = («!,..., ctrf) and a number, — a^-i-i, called the offset. Now, a hyperplane h consists of all 
points X G E'^ such that 

(x, a) + ad+i = 0, (5-a) 

that is, the scalar product of x and a equals the offset. Notice that the hyperplane does not 
change if we multiply the normal vector and the offset by some nonzero number. We define h* as 
the point whose homogeneous coordinates are (cti, . . . , a^; ctd+i). Geometrically speaking, h* lies 
on the line through the origin defined by a, and the distance of h* from the origin is the inverse 
of the distance between h and the origin. This can easily verified after observing that is 
the distance between h and the origin, provided a has unit length. Note also that the origin 
lies between h and h* (see Figure 5-II). Conversely, for a point p with homogeneous coordinates 
(tti, 7r2, . . . , TT^; TTrf+i) we let p* be the hyperplane with normal vector (tti, 7r2, . . . , tt^) and offset 




Figure 5-II: Mapping a line to a point and vice versa. 

It is straightforward to show that this transformation preserves incidences; that is, p G /i if and 
only if h* G p*. Indeed, it is a triviality when one remembers what p ^ h means algebraically, 
namely, that 

It is equally easy to prove that this mapping preserves the relative order between a point and a 
hyperplane. To describe what exactly we mean by this define 

= {x\(x, a) -\- ctd+i > 0} and h~ = {x\(x, a) -\- a^+i < 0}, 
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and call those the positive and negative sides or half-spaces of h. By order preservation we mean 
that p if and only if h* G p*~^ ■ Here, a warning is appropriate to avoid future confusion. 

If we multiply the normal vector and the offset of a hyperplane h hy —1 we do not change the 
hyperplane but we do change the sides of h; what was previously its positive side is now its negative 
one and the other way round. We will take advantage of this curiosity by encoding the positive 
and negative sides into the hyperplane's specification. Note that geometrically, the normal vector 
of a hyperplane h points to its positive side. 

The primitive operation that we wish to tackle in this section is to decide on which side of a 
hyperplane hi^ the intersection of d other hyperplanes hi^, . . . , ^id-i ^i^^- -^y °^ SoS, the 

absence of any kind of degeneracies can be assumed; so hi^ (e) through hi^_^ (e) intersect in a unique 
point which does not lie on hi^[e). By Cramer's rule, the intersection point p = (tti, 7r2, . . . , tt^) of 
d hyperplanes is given by the coordinates 



det A 



TT,: 



det A„ 



where is the matrix 



and Ad^i is the same matrix after replacing the z-th column from the left by the vector 



/ 



\ -a»d-i,d+i / 

Point p lies in the positive half-space of if and only if 

Provided that det A^ is positive, this is equivalent to 

det Ad,ia,^^i + det Ad,2ajd,2 H h det Ad,da,^^d + det Ada,^^d+i > 0. 

In case of a negative det A^ the above statement is valid after reversing the direction of the 
inequality. Consequently, p €^ hj^ if and only if 



det Ad+i ■ det A^ > 0. 



This can be seen by developing 



det Ad+i = det 



««o,2 
ttjl ,2 



a 



«d-l ,2 
««d,2 



^io :d 



Ciio,d+l 
Ctil ,d+l 



a; 



Id- 



-i4 «• 



4+1 



using the last row. Now, we can use this to write a procedure that decides on which side of a 
hyperplane J other hyperplanes intersect. It uses SoS, as described in Section 3.3. 
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Predicate 4 (OnPositiveSide) Let hi^^hi^, . . . ,hi^ be J + 1 hyperplanes in J dimensions, given 
as in (5-a), and with distinct indices < Zq, «i, • • • , < — 1. The following function written in 
pseudocode returns true if the intersection 0^=0 ^u(^) ^^^^ positive half-space of hi^[e) and 

false if it lies in the negative one. 



function OnPositiveSidedihia, ■ ■ ■ ,hi^_^]hi^) returns Boolean 

loral i' i' d' i" i" i" <i" d" 

begin 

Sortd{{io, • • • , irf-i), (i!o5 • • • , s') 
Sortd+i{{io, • • • , irf-i, id), (i!o, • • • , id-i,id)^ s") 



d' = SignDetAa 



if odd (s') then d' 



d" <— SignDetAd+i 



if odd (s") then rf" 
return (d' = d") 
end 



'd-i' 

a ' 



"«J,,2(e) • • 








".■J,',2(e) 
"«;',2(e) 


• • • 

• • • 


".■^i,2(e) 
"«;;,2(e) 


• • • cij" (i(e 

• • • 



Q!^^',d+i(e) / 



5.3 Nonvertical Hyperplanes 

In many applications we know that all hyperplanes we have to deal with are nonvertical; that 
is, they intersect the J-th coordinate axis in a unique point. Examples are Voronoi diagrams 
or, more generally, power diagrams for arbitrary order and weighted Voronoi diagrams (see for 
instance [Ed87] and [AI86]). It is beyond the scope of this paper to describe how the data for those 
problems are used to generate hyperplanes — it will be enough to know that they are obtained 
via geometric transforms which do not create vertical hyperplanes. 

A nonvertical hyperplane h in J-dimensions can be specified by a relation of the form 

aixi + a2X2 H h ad-iXd-i + Xd + ad = 0. (5-b) 

The advantage of describing a hyperplane using this form rather than the one in Section 5.2 is 
that it takes only d parameters rather than d-\-l. This will lead to some savings when it comes to 
computing signs of determinants (compare for instance det A4 in Table 6 and det A4 in Table 5). 
Since every hyperplane h is now nonvertical, we can uniquely define what we mean when we say 
that a point lies (vertically) above or below h. Define 

/i+ = {x = (xi, . . . , Xd)\aiXi H h ctd-i + Xd + ad> 0}, 

and let h~ = E"^ — h — . A point p is said to lie above h if p €^ and below h if p €^ h~ . 
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The primitive operation that we consider in this section decides whether the intersection of d 
hyperplanes hi^, . . . ,hi^_^ Res above or below hyperplane hi^. The use of SoS as in Section 3.2 
allows us to assume that indeed hi^ through hi^_^ intersect in a unique point that does not lie on 
hi^. A decision procedure based on comparing the signs of the two determinants can be derived 
from the procedure given in Section 5.2. We just replace all ai^di^) by 1 and exchange the last two 
columns of the second matrix in function OnPositiveSide d- This leads to the following predicate. 

Predicate 5 (^Ahove) Let hi^^ hi^, . . . , hi^ he d -\- 1 nonvertical hyperplanes in E'^ specified as in 
(5-b), and with pairwise different indices < Zq, «i, . . . , < — f . The following predicate returns 
true if the point of intersection 0^=0 ^u(^) ^^^^ above hi^[e)^ and false if it lies below hi^[e). 



function Aboved(hia , . . . , hi^_^ ; /ij-^) returns Boolean 

loral ■}' i' d' i" i" i" <i" d" 

begin 

Sortd{{io, • • • , id-i), {i'o, i'd-i), s') 
Sortd+i{{io, • • • , id-i,id), {i'o, • • • , i'd-i^i'd), s") 

d' = SignDetAd \ ' ■ . \ \ 

if odd (s') then d' < d' 



d" <— SignDetAd+i 



a^'' d_i(e) a!j''d(e) 1 \ 



cij" lie I 



a^" d_i(e) a!j''d(e) 1 / 



if odd (s") then d" 
return (d' ^ d") 
end 



d" 



5.4 In-sphere Test 

In J dimensions any J + I affinely independent points (i.e., points that do not lie in a common 
hyperplane) define a unique sphere that goes through the J + I points. For example, in two 
dimensions there is a unique circle through any three noncoUinear points. Given J + 2 points 
po,Pi, • • • ,Pd+i the problem we address in this section is how we can determine whether pd+i lies 
inside or outside the sphere specified by the first J + I points, assuming this sphere is unique. 
Such a test is useful for constructing Voronoi diagrams (as shown in [GS85] for d = 2) and other 
problems where circles and spheres play a role. 

An elegant solution to this problem can be given using a transform that lifts a sphere in J di- 
mensions to J + I dimensions where it is represented by a hyperplane. This transformation can 
be traced back in the literature to [Se82] and has since been used throughout the computational 
geometry literature (see [GS85], [PS85], and [Ed87]). For the case of circles in the plane we explain 
this transformation in detail and finally phrase the predicate for general dimensions. 

Let U : X3 = xf A X2 he the paraboloid of revolution whose symmetry axis is the Xs-axis, and let 



c : {xi- 7i)2 + {x2 - 72)^ = 73 
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be a circle in the a;ia;2-plane. Note that (71,72) is the center of the circle and 73 is its radius. The 
lifting map transforms c to the plane c* in three dimensions given by the equation 



2^ 



c* : X3 = 271X1 + 272X2 - (7i +72-73 

The quick reader will already have verified that the vertical projection of U He*, which is an ellipse 
in three dimensions, onto the a;ia;2-plane is equal to the original circle c. A point p = (7ri,7r2) lies 
inside c if and only if its vertical projection onto U lies below c*. This insight gives us some hope 
that, in fact, the problem can be bent such that Predicate 5 from the previous section is applicable. 
Before we continue our exploration in this direction, let us understand how the original statement 
of the problem and the lifting map are connected. Recall that there are four original points, 
which we call po, Pi, P2, and ^3. The first three determine the circle c and therefore the plane 
c*. Moreover, if we project them vertically onto [/, then c* is the plane through these points on 
the paraboloid U. The question is now whether p'^ = (vr3^i, 7r3^2, 7r| + vr|2), which is the vertical 
projection of ^3 onto [/, lies below c* (in which case ^3 lies inside c) or above c* (then ^3 lies 
outside c). By the use of SoS, we can assume that the four points are in general position. 

This problem can be mapped to the plane problem of the previous section if we use a dual 
transform. This transform replaces each point on U by the unique plane whose intersection with 
U is this point. If p'^ = (tti, 7r2, irf + ^2)^ then the formula for this dual plane is 

p* : X3 = 27ria;i + 27r2a;2 — (vr^ + tt^). 

We see that this is indeed the lifting map applied to point p = (7ri,7r2) in the a;ia;2-plane. This 
duality transform preserves incidences and above-below order in a way similar to the duality 
transform described in Section 5.2. This leaves us with the following correspondence between the 
original point-circle problem and the derived plane-point problem: Point ^3 lies inside c (the circle 
through points po, Pi, and P2) if and only if the intersection point of the planes p^^ pj, and P2 lies 
below P3. The statement is also valid if we replace "inside c" by "outside c" and "below ^3" by 
"above ^3" . 

We leave the generalization of this two-dimensional exercise to three and higher dimensions to 
the curious reader. In any case. Predicate 5 can now be used to implement Predicate 6 which 
formalizes the in-sphere test in d dimensions. If we apply Predicate 5 directly, we will find ourselves 
computing the sign of determinants of the form 

/ 27r,„,i ••• 27r,„,, -(^r^ + . . . + vr^ J -1 \ 



det 



The sign does not change if we divide the entries in the left d columns by 2. Similarly, we can 
remove the minus signs in the last two columns without changing the sign of the determinant. 
However, there remains one problem with determinants of the above type, and this is that the 
values in the [d + l)-st column from the left depend on the values in the left d columns. In 
particular, with SoS, the e-expressions of the point coordinates appear in mixed products in the 
(d + l)-st column. This turns out to be a real pain when we implement SoS for this type of 
determinants. A cheap trick that handles this problem is not to perturb the original points but 
rather to perturb the vertical projections onto the paraboloid in J + 1 dimensions. In effect, this 
means that we introduce 

d 

7r.„d+i = J2<,. for 0<A<J+1, (5-c) 
i/=i 
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and then perturb the points (vri;^,i, • • • , TTip^^^, TTi^.d+i)- Because the perturbation of the [d + l)-st 
coefficient does not depend on the first d coefficients, this imphes that the points are perturbed 
away from the paraboloid U. On the other hand, if the perturbation is small enough we are still 
close enough to the original situation. 



Predicate 6 (InSphere) Let pi^^pi^^ . . . ,Pi^j^^ be J + 2 points in J dimensions with pairwise dif- 
ferent indices in the range from through n — 1. The program below returns true if the perturbed 
image of Pi^_^_^ lies inside the sphere through the perturbed images of the first d -\- 1 points, and 
returns false if it lies outside. 



function InSphere^ (pi^^, . . . ,pi/,pi^_^^) returns Boolean 

loral ■}' i' d' i" i" i" <i" d" 

begin 

SortdJ,i{{iQ, id), {i'o, . . . , i'^), s') 
Sortd+2{{io, ■■■,id, id+i), {i'o, i'd, i'd+i), s") 

d' = SignDetAd+i '■ 



vr,-' Je) 1 



if odd (s') then d' < d' 

Set TTi^^d+i as in (5-c). 



d" ^ SignDetAd+2 



if odd (s") then d" < d" 

return (d' = d") 



end 



Note that the rightmost column of the first matrix in the above program should really consist of 
— I's. To stress the similarity with predicate Abovcd+i in the previous section we replaced the — I's 
by -|-l's and thus changed the sign of d' . This effect is compensated by the fact that we want to 
return true where function Abovcd+i returns false. 



6 Remarks and Discussion 



The main contribution of this paper is the introduction of a general technique that can be used 
to deal with degenerate input for geometric programs. The main purpose of this paper is to 
demonstrate that this technique (which we call SoS, the Simulation of Simplicity) is immanently 
practical, despite its high-powered appearance. Indeed, the authors believe that SoS will become a 
standard tool for implementing geometric algorithms. A pragmatic consequence of this technique 
is that authors of geometric algorithms can now be more confident about the implementability of 
their algorithms even in the presence of any conceivable degeneracies, provided SoS is applicable 
to their algorithms. 
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This raises the question of determining the hmitations of SoS — what are the properties of an 
algorithm that allows us to use SoS when we implement it? One important feature of algorithms 
that are amenable to SoS is that their algebraic computations are of constant depth. The deeper 
the algebraic computation the more complicated is the polynomial (or, in general, the function) 
in e generated by SoS, and the less tractable is its evaluation. Another limitation of SoS is the 
necessity of absolute precision in the evaluation of algebraic formulas. As long as square roots 
can be eliminated by squaring the equation and similar techniques can be used to remove other 
irrational functions this is not a problem, but there are cases where it is not that easy. Typical 
examples for such problem cases are algorithms for shortest path problems in a geometric setting. 
Take for instance two piecewise linear paths in the Euclidean plane. The length of each path is 
the sum of square roots of integers (assuming the endpoint coordinates are integers). Deciding 
which one of the two paths is shorter is a difficult question unless the number of square roots is 
very small. On the other hand, deciding which one of two paths is shorter is not exactly the kind 
of problems that SoS was invented for. 

Another problem is that algorithms employing SoS produce results for the perturbed set of objects 
rather than for the original ones. In certain settings, such as in computer graphics, this fact can 
often be ignored. However, when "unperturbed" results are needed, some postprocessing has to be 
performed. This paper does not deal with this issue and further work has to be done. Nevertheless, 
in most of the applications mentioned in this paper the postprocessing step is more or less trivial: 

- In the point-in-polygon problem one can simply add a test whether or not the query 
point lies on a boundary edge. 

- In the case of Voronoi diagrams or arrangements of hyperplanes, we identify and elim- 
inate zero-length edges or higher-dimensional faces of zero measure. 

- In the convex-hull setting, it is possible to undo the perturbation simply by merging 
adjacent faces if necessary; for example, in two dimensions, adjacent edges that lie on 
a common line, and in three dimensions, adjacent triangles contained in a common 
plane. 

It is rather difficult, however, to use SoS or any other perturbation scheme for finding all data 
points on the boundary of the convex hull. This is because the perturbation may decide that 
a point is inside the hull if it lies on a boundary edge or face. In this case the point would be 
prematurely discarded. We refer to [Ya87] for a more extensive discussion of the limitations of 
symbolic methods aimed at resolving robustness problems in geometric algorithms. 

In order to increase the credibility of our claim that SoS is indeed a practical programming tool, the 
second author compiled a prototype version of a SoS library [Mu88] and implemented the three- 
dimensional edge-skeleton algorithm of [Ed86]. We believe it is fair to say that this algorithm is 
an extraordinary challenge for someone who wants to do it without SoS. From run-time profiles of 
this program we learned that most of the computing time was spent on multiplying long integers 
in order to compute signs of determinants. The speed-up that we got in our implementation 
from replacing long-integer by normal (built-in) integer arithmetic was a factor somewhere around 
10. Of course for the normal integer arithmetic to work we severely restricted the range of the 
coordinates that were used. In any case, this makes it clear where future work has to go if we 
want to produce programs that are reliable and which are as fast as software that uses fioating- 
point arithmetic and is therefore inherently unreliable. The most promising way to eliminate this 
overhead factor seems to be the design of a special piece of hardware that computes the sign of 
determinants for integer matrices. Such effort seems justified by the versatility of determinants 
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demonstrated in Section 5. We would like to mention, though, that even without the availability 
of such specialized hardware we believe that SoS is of practical value in implementing geometric 
algorithms. Aside from the obvious savings in time and effort for the programmer, it seems to us 
that the use of SoS is currently the only hope to produce geometric software that is in any sense 
reliable. 

We end this section by pointing out a new direction for further research — it is the systematic 
study of primitive operations used and needed for geometric algorithms. If one undertook the 
venture of building a library of primitives for geometric algorithms, besides computing signs of 
determinants, what other operations would have to be in the collection? Is it even clear that 
computing the sign of a determinant is such an indispensable operation or are there less expensive 
ways to determine the orientation of J + 1 points in d dimensions? 
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Appendix 

In this appendix, we give the relevant subdeterminants, sorted in sequence of decreasing signif- 
icance, needed for computing the signs of det A4(e), detA2(e), det A3(e), and detA4(e). Each 
sequence is given in a table that also shows the corresponding e-product Ct and the size kf of the 
matrix Mf^'' [M^"^) associated with the it + l)-st significant term in the e-polynomial det Aj:)(e) 
(det Aj:)(e)). The third column of each table shows Uj, the vector that encodes the subdeterminant 
of depth t. Recall that this vector was used to produce the proper sequences of subdeterminants 
by successive calls of procedure Next-V. 
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kt ■ kt 




det 2 







2 • 2 


[2, 2; 2] 




e{) 
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1 • 1 


[1,2; 2] 


+ det(l) = +1 


e(^,l) 



Table -i: The 2 relevant terms of det A2(e). 
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Table -ii: The 5 relevant terms of det A3(e). 
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Table -iii: The 15 relevant terms of det A4(e). 
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Table -iv: The 50 relevant terms of det A4(e) (to be continued). 
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Table 6 continued. 
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Table 6 continued. 



